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which  can  be  efficiently  Implemented  directly  in  floating-point  arithmetic  with  no 
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CHAPTER  I 

Introduction 


Let  a and  k be  positive  integers  and  let  p be  a prime.  Let  J{x)  • xk-ak,ixk~} a0  be 

a primitive  polynomial  in  Zp[x],  In  other  words,  the  residue  classes  modulo  J{x)  of  the 
polynomial  ring  Z^x]  form  a finite  field  in  which  x Is  a generator  of  the  multiplicative 
group.  Consider  the  Jk‘h-order  linear  homogeneous  recurrence 

JW«A-iWi  + ‘ • +«o7nmod^a  01) 

for  n - 0,  1,2,  ...  and  initial  values  (y^  ....  jA_,)  * (0,  ....  0)  mod  p.  The  sequence  of 
fractions  <y,]pa>  is  a candidate  for  a pseudo-random  sequence.  If  a - I then  the  length 
of  the  period  of  <yn>  is  pk- 1,  and  if  a>  '■  the  length  of  the  period  is  pa~c(pk- 1),  where  c Is 
easy  to  compute  and  is  often  equal  t - Recurrences  of  type  (1.1)  with  a - * - I and  large 
p,  or  with  k - I,  p - 2,  and  moderately  '*rge  o,  are  the  basis  for  some  of  the  most 
acceptable  and  widely  used  computer  random  number  generators  [Knuth69X  Such 
generators  have  |a0| > I and  require  that  a multiplication  modulo  pa  be  performed. 
Multiplication  can  be  replaced  by  addition  and  subtraction  if  I;  of  course  then  k must 
be  greater  than  I.  Experimental  evidence  has  been  accumulating  to  the  effect  that 
recurrences  of  type  (l.l)  with  p - 2 and  moderately  large  a and  k are  quite  successful.  (See 
[Knuth69]  p.  464  and  [Brent73]  pp.  163-164;  also  [Green59]  and  [Franklin64X>  However, 
theoretical  Justification  for  such  success  has  been  lacking.  We  will  show  that  recurrences  of 
type  (l.l)  are  indeed  excellent  random  number  generators  by  showing  that  for  tsk  the 
Muples  of  consecutive  residues  become  equidistributed  in  t dimensions  in  the  limit  as 
a -*  m,  subject  only  to  a much  weaker  condition  on  the  distribution  of  the  residues. 

In  the  remainder  of  this  chapter  we  shall  discuss  known  results  on  the  length  of  the  period 
of  sequences  of  type  (l.l).  Chapter  2 defines  discrepancy,  a means  of  measuring 
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equidistribution,  and  presents  a formula  of  Harald  Niederreiter  which  expresses  the 
discrepancy  in  terms  of  an  exponential  sum.  In  Chapter  3,  exponential  sums  are  used  to 
reduce  the  question  of  equidistribution  of  sequences  of  type  (l.l)  to  a much  weaker 
distribution  criterion,  and  in  Chapter  4 the  sequences  are  analyzed  with  respect  to  this 
criterion.  It  is  believed  that  the  analysis  of  Chapters  3 and  4 is  new.  Chapter  5 considers 
implementation  details  and  gives  the  results  of  empirical  tests  comparing  higher  order 
linear  congruential  generators  of  type  (i.i)  with  standard  linear  congruential  generators 

yn*i~ayn*bmod  Pa- 

The  simplest  example  of  a sequence  satisfying  (l.l)  with  Jk>l  is  the  Fibonacci  sequence 
with  pm  2 and  initial  conditions  (y0,  ys)  * (0,  i).  The  recurrence  is  Fn+2  « Fn+,  + Fn, 
corresponding  to  the  primitive  polynomial  x2-x-l  - x2  + x+  I in  Z2txl  The  period  is 
3.2a'‘  and  it  is  known  that  the  sequence  of  ordered  pairs  (2'aF n 2'aFn+1)  becomes 
evenly  distributed  mod  i as  a increases.  (See  [Marsaglia72).)  However,  the  Fibonacci 
sequence  is  not  a suitable  random  number  generator  because  successive  triples  are  very 
poorly  distributed  in  three  dimensions.  To  achieve  satisfactory  performance  we  must 
consider  recurrences  of  higher  degree. 

When  considering  such  recurrences  it  is  helpful  to  know  some  facts  about  the  length  of  the 
period  and  some  relationships  between  sequences  satisfying  the  same  linear  congruence,  but 
with  different  initial  conditions.  The  papers  [Ward3IJ,  [Ward33),  and  [Hal!38a]  present 
accounts  of  the  theory  for  general  linear  recurrences.  The  length  of  the  period  of 
recurrence  (l.l)  for  a - I can  be  easily  established  using  an  idea  from  [Hall38a]. 

Lemma  l.l.  If  fix)  is  primitive  in  Zp[x]  then  the  period  of  (l.l)  for  a » I is  pk-  I. 

Proof.  Corresponding  to  the  Jk-tuple  (y^  . . . , yA. ,),  associate  the  polynomial 

•••  0-2) 
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Then  xK(x)  - *,**  + (*!-«*., J0>**',+  *"  +(?a-i aiJo>*  - + (J2-«*-iJi>*A'2 

+ •••  + (jA a^,)  (modulo /(x)),  and  this  is  the  polynomial  associated  with  (y, 

yk).  Thus  iterating  the  recurrence  corresponds  to  multiplying  the  associated  polynomial  by 
x.  Since  f{x)  is  primitive,  x is  a generator  of  the  multiplicative  group  of  the  finite  field 
Zp[x]l(J{x))  and  has  period  pk-  I.  Thus  recurrence  (1.1)  has  period  pk~  I.  I 

We  will  use  generating  functions  and  congruences  to  a double  modulus  to  analyze  the 
period  and  certain  other  properties  of  integer  sequences  satisfying  (1.1)  when  a>  I.  The 
notion  of  congruence  to  a double  modulus  is  an  extension  of  the  usual  notion  of 
congruence.  If  J{x),  a(x),  and  b(x)  are  polynomials  with  integer  coefficients,  and  if  there 
exist  polynomials  u(x)  and  v(x)  with  integer  coefficients  such  that 
a(x)  m b(x)+f{x)u(x)  + mv(x),  then  we  will  write  a(x)  m b(x)  (modulo /(x)  and  m),  or  also 
a(x)  s b(x)  (modd  m,  /(*)).  The  following  four  lemmas  are  essentially  exercise  3.2.2-11  of 
[Knuth69l 

Lemma  1.2.  Assume  that/(0)  is  relatively  prime  to  p and  that  pa>  2.  If 

xr  m I (modd  pa, fix))  and  xT  a I (modd  pa*\/(x)) , 

then 

x*  k I (modd  pa  * \f[x))  and  x**  * I (modd  pa  * 2,  f{x)) . 


Proof.  By  definition  of  congruence  to  a double  modulus  there  exist  polynomials  u(x),  v(x) 
such  that  xT  - I +f(x)u(x)+ pav(x).  In  addition  o(x)  a 0 (modd  p,  f{x)),  or  else 
xT  - I + /(x)u(x)  + pa(J{x)uz{x)  + pv2[x))  - I +/(xX«(x)+ pau2(x))+  pa  * *»a(x)  ■ I (modd 
P°  * /(*)).  contradicting  the  assumption  about  xT  (modd  pa  * J[x)).  If  we  raise  xr  to 

the  pth  power,  the  binomial  theorem  makes  it  clear  that 
l + pa*lv(x)+p2a*lvz(x)(,p-\)f2  plus  other  terms  which  are  congruent-  to  zero 
modulo /(x)  and  pa**.  Since  pa> 2 we  have  pka*i(p-l)l2  a 0 modulo  pa*2,  and  thus 
x^  ■ I + pa  * Mx)  (modd  pa  * *,  J{x)).  Suppose  that  pa  * Vx)  ■ 0 (modd  pa*z,  /(x)). 


I 
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Then  there  exist  polynomials  a(x),  b(x)  such  that  pa  * *v(x)  - a(x)/(x)  + pa  * zb{x),  which 
implies  that  pa  * 1(v{x)-pb{x))  - a(x)J{x).  Since  /( 0)  is  relatively  prime  to  p we  can  apply 
Gauss’s  lemma  about  the  g.c.d.  of  the  coefficients  of  the  product  a(x)J[x)  and  deduce  that 
pa  * 1 divides  a{x).  This  means  that  v(x)  - yf (a + 1 ^a{x)f[x)  + pb(x)  s 0 (modd  p,f{x)),  which 
is  a contradiction.  Thus  pa  * V*)  * 0 (modd  pa*z,  fix))  and  therefore  * I (modd 
pa'Z,fix)).  I 

Using  generating  functions  we  will  next  derive  a relationship  between  period  lengths  and 
the  powers  T for  which  xr  s I (modd  m.fix)).  Let/(x)  - l-o,x-  • • • -akxk,  and  let  G(x)  - 
Mfix)  - A0+  >4,x  + A#*  • • • . Denote  by  T(m)  the  length  of  the  period  of  <An  mod  m>. 

Lemma  1.3.  T(m)  is  the  least  positive  integer  T such  that  xT  s 1 (modd  m.fix)). 

Proof.  Since  T (m)  is  the  period  of  <A„  mod  m>  we  have  G(x)-xr^m^G(x)  e 
A0+Atx  + A^xz  + ••  ^T(m)-txT(m)1  (modulo  m),  which  implies  l-xT(m)  = 
fix)(A0+  Ayx  + •••  +>*T(m)-ixT(m)  *)  (modu,°  m)>  80  l-xT(m)  s 0 (modd  m,  fix)).  This 
shows  that  T(m)fcT,  since  T was  the  smallest  positive  integer  such  that  I -xT  e 0 (modd  m, 
fix)).  Conversely,  taking  xT  - I e 0 (modd  m,fix))  and  multiplying  by  G(x)  - 1 I fix)  gives 
xTC(x)-C(x)  = 0 (modulo  m).  Equating  co-fficlents  of  x gives  An-An^T  £ 0 (modulo  m) 
for  all  n*0.  Thus  T(m)sT.  I 

We  now  restrict  our  attention  to  prime  power  moduli  m - pa  and  show  that  when  a Is 
large  enough,  increasing  a by  I multiplies  not  only  the  modulus  but  also  the  period  by  p. 

Lemma  1.4.  If  pa  > 2 and  T(pa)  * T(pa  * *)  then  T (pa  * k)  - pk7{pa). 

Proof.  It  suffices  to  show  that  T{pa)  * T(pa  * ’)  implies  that  T{pa  + ')  - pT(pa)  * T(pa  * z). 
From  what  we  have  already  shown  about  the  period,  we  know  that  r(pa  * 2)  * pT(pa)  and 
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that  T (pa  * ')  divides  pT(pa)  but  does  not  divide  T (pa).  Let  T (pa)  - p*q  with  p and  q 
relatively  prime.  Then  since  T(pa  * *)  divides  p**lq  but  does  not  divide  peq,  it  must  be 
that  r(pa*  ’)  <«  p**xd  where  d divides  q.  Since  p**xd  is  a period  modulo  pa*  *,  it  is 
certainly  also  a period  modulo  pa.  The  smallest  period  divides  ait  other  periods,  so  pmq 
divides  p**ld.  Thus  q divides  d,  q - d and  t(pa  4 *) » pi{pa).  I 

We  have  determined  the  period  of  <An  mod  pa>  where  An  is  the  coefficient  of  x"  in  the 
generating  function  G(x)  - Mf{x)  and  f(x)  - I -a,*-  ••  • -aAx*.  This  sequence  <An  mod 
pa>  satisfies  (l.l),  and  we  would  like  to  know  something  about  the  set  of  all  sequences 
satisfying  (1.1).  If  <?,(*)-  J0+J|*  + where  yn. k - a*.,yntJk.,  + •••  +a(pa  for  all  naO, 
it  is  easy  to  see  that  f{x)G(x ) is  the  polynomial  gix) m y0  + (yl-aty0)x+  ■ • • 

+ (?*-! ~ai?A-2_a2?A-3~ flA-lJo)*A  hence  glxMx)  ,s  the  generating  function  for 

the  sequence  with  initial  values  (y0,  yv  ....  The  next  lemma  shows  that  the  period 

of  <ya  mod  m>  is  the  same  as  the  period  of  <A„  mod  m>,  in  the  cases  of  interest  to  us. 

Lemma  1.5.  Let  m - pa.  If  f(x)  and  g(x)  are  relatively  prime  modulo  p then  the  period  of 
<yn  mod  m>  equals  the  period  of  <An  mod  m>.  In  particular,  this  holds  when  J{x)  is 
irreducible  and  (Jq,  . . . , yk. ,)  * (0 0)  (modulo  p). 

Proof.  Assume  that  T is  the  period  of  <yn  mod  m>.  Then  g(xH I -xT)  h 0 (modd  m, /(*)). 
Because /(x)  and  g(x)  are  relatively  prime  modulo  p,  we  can  apply  Hensel’s  lemma  and  find 
polynomials  a(x)  and  <Kx)  such  that  a(x)f{x)+b(x)g{x)  s I mod  m.  Multiplying 
g(xKl  -xT)  b 0 (modd  m,  /(x))  by  <Kx)  gives  I -xT  s 0 (modd  m,  /(x)),  hence  the  period  of 
<yn  mod  w>  is  no  shorter  than  the  period  of  <An  mod  m>.  On  the  other  hand  it  cannot  be 
longer,  since  <yn  mod  m>  is  a linear  combination  of  sequences  <An+j  mod  m>  for  various/ 
If  f(x)  is  irreducible  in  Zj[xJ  then  f(x)  is  relatively  prime  to  every  nonzero  polynomial  of 
lower  degree,  so  the  period  of  <yn  mod  m>  will  be  the  same  as  the  period  of  <An  mod  m> 
unless  g<x)  s 0 mod  p.  This  can  only  happen  if  all  of  the  initial  values  (yp . . . £a.j)  are 
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divisible  by  p (again  because  of  Gauss’s  lemma).  I 


The  results  which  have  been  derived  so  far  about  the  period  may  also  be  derived  by 
formulating  the  recurrence  relation  (1.1)  in  terms  of  a matrix-vector  product.  This 
alternative  formulation  is  important  because  actual  computation  may  be  easier  with  the 
matrix  than  with  the  polynomials.  Let 


0 1 0 - - 

0 0 1 ••• 


0 0 0 ••• 
Uoaia2  ••• 


(1.3) 


where  the  are  the  same  coefficients  as  those  in  (1.1).  Define 

- (yn,  yntV  . . to  be  the  *-tuple  of  consecutive  terms  of  (1.1)  beginning  with 

the  nth  term.  It  is  easy  to  see  that  s Ay0  and  hence  by  induction  that  yn  s A nyQ  mod 
pa.  By  definition  of  the  period  T - T(p°)  we  have  jn  s Atj0  mod  pa,  so 

At  ■ I*paB  (1.4) 

for  some  matrix  B.  If  p does  not  divide  B then  T is  not  the  period  mod  pa  * V However, 
A^  ■ (1  + paB)p  - 1 + pa  * ’B  + P'  ^^-•pZaB2  plus  other  terms  divisible  by  pa*2.  If 
pa  >2  then  we  see  that  A ^ s 1 + pa  * *B  (modulo  pQ  * 2),  and  we  have  rederived  the  fact 
that  pr  is  the  period  modulo  pa*x  but  not  mochi'o  pa*z.  We  can  also  note  that 
jn+T/p-Jn  s AT/pjn-j(n  E P**  'fyn  (modulo  pa ).  This  says  that  as  the  modulus  increases 
the  difference  between  vectors  a period  apart  is  merely  multiplied  by  p.  In  the  binary  case 
p - 2 this  means  that  the  exdusive-or  of  yntj/2  and  jn  gives  the  same  pattern  of  bits,  just 
shifted  one  place  left  as  the  modulus  increases. 

If  we  consider  A^  a}  for  j ■ pa  we  see  that 

A^'T<a>-(UpaB/,-I*papaB^a*^f2aB2+  - - » I + p^B  mod  p3a‘l  (1.5) 


Introduction 


When  Eq.  <1.4)  is  viewed  as  defining  B as  a function  of  a then  Eq.  (1.5)  says  that 
B(a)  s B(2a)  (modulo  pa‘x).  In  other  words,  as  a increases  B converges  in  the  p- adic 
sense. 
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CHAPTER  2 

Discrepancy  and  Exponential  Sums 

In  order  to  make  meaningful  statements  about  how  well  a sequence  Is  distributed  It  Is 
necessary  to  have  a measure  of  equidistribution.  Discrepancy  is  such  a measure.  Let  I be  a 
unit  interval,  and  for  intervals  J = I let  A(J,Af)  be  the  number  of  n,  Osn<N,  with  yn€ J. 
Then  the  discrepancy  of  the  points  j#  ....  jN.x  is  defined  as 

• • • ’In-*)  ” I ^ (2.1) 

where.  A(J)  is  the  measure  of  J.  Thus  the  discrepancy  measures  the  maximum  difference 
between  the  actual  fraction  of  hits  in  an  uterva!  and  the  expected  fraction  of  hits.  It  is 
easy  to  see  that  OiDNi  I and  that  the  sequences  produced  by  a good  random  number 
generator  should  have  small  discrepancy.  Definition  2.1  can  be  extended  naturally  to 
define  discrepancy  for  sequences  of  points  %n  lying  in  a multidimensional  unit  interval  I, 
and  we  shall  often  use  the  extended  definition. 

Harald  Niederreiter  has  developed  an  inequality  relating  the  discrepancy  to  certain 
exponential  sums  [Niederreiter 78],  This  inequality  is  important  because  it  bounds  the 
discrepancy  in  terms  of  functions  having  nice  mathematical  properties.  The  properties  will 
be  exploited  when  analyzing  the  discrepancy  of  linear  congruential  sequences  (l.l). 

Niederreiter’s  inequality  Is  easier  to  state  if  some  notation  is  introduced.  For  integers  m and 
A define 

I if  A * 0 (mod  m) 

r(A,  m)  - 

m sin  n||A/fn|)  if  A a 0 (mod  m) 
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where  ||/ 1|  is  the  distance  from  t to  the  nearest  Integer.  For  lattice  points 4 - (A,,  ....  hB)  * 
Za  we  write 


r(im)-  FI 

The  summation  symbol  2)  designate  a sum  over  the  complete  system  of 

h mod  m 

numerically  least  residues  modulo  m,  consisting  of  all  integers  A with  -m/2< Asm/2.  The 
* 


summation  symbol 


E 

h mod  m 


refers  to  the  same  sum  but  with 
* 


A 


0 deleted  from  the  range  of 


summation.  The  symbols  2 and  2 wil1  refer  to  analogous  sums  over  the 

h mod  m h mod  m 

complete  system  of  representatives  of  Za/(mZ)a  consisting  of  all  h * (A, A#)  c Z*  with 

-m/2< A^Sm/2  for  I <jis,  possibly  omitting  A - (0,  ....  0) . The  notation  jc*y  represents 
the  standard  inner  product  of  two  vectors.  The  function  e(t)  is  defined  for  real  values  t as 
e(t)  - . 


Lemma  2.1  (Nlederreiter’s  lemma).  Let  Jq,  be  N lattice  points  in  Z*.  Then  for 

any  integer  m^2  the  discrepancy  DN  of  the  points  m"1^  ....  m~'%N_ , satisfies 


Dn< 


s 

m 


* 

♦ S 

h mod  m 


I 

K 4,  m) 


Proof.  For  4 - (A, kt)  * Z*.  let  A(t;A/)  - A(A, ka;N)  be  the  number  of  n, 

Os  ns  N-  1,  such  that  = 4 (modulo  m),  and  let  ek  be  the  characteristic  function  of  the 

coset  4+ (mZ)*  of  Za/(mZ)a.  Then  for  g - (x ,,  ....  *,) « Z*  we  have 

3 

t{hpCj-k j)lm) ) ■ —B  2 •(AKs - 4 )/»»)  • 

~ W Jm\  hj  mod  m m h mod  m 

Therefore, 


N- 1 N- 1 . N- J 

A (4,  A/)  ■ 2 — 2 £ — 2 *< “ A*A/«)  2 • 

n*0  — W n*0  /»  mod  m m h mod  m n-0 

and  so 
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Mk;N)- —t  • ± 2 <<-4*|/«)2  <2.2) 

W Vt  h mod  m a*  0 

Now  let  J - la,,  0,)  x ...  x la,,  (3,)  be  an  arbitrary  half-open  subinterval  of  [0,  I)4. 
For  each  J,  Isjss,  we  choose  the  largest  closed  subinterval  of  C ajfij)  of  the  form  [ujm, 
v Jut ] with  integers  Uji  Vj.  The  case  where  for  some  J no  such  subinterval  of  t aj,  0j)  exists 
can  be  dealt  with  easily,  since  we  have  then  A(J,A0  » 0 and  0 j-Oj<  ^ . hence 


|^^-V(J)|  -V(J)<  1. 


In  the  remaining  case,  the  integers  u,,  . . . , . ...  vt  are  well  defined,  and  we  obtain 


A(J;A0-AJV(J).  2 (a(4;/V)-  -^)  4 -^(oj-Uj+I)  •••  (p,-u,+  1)  -AfV(J) 

ktujikjivj  m m 


■ — 2 ( Tj  4-4-4/w))(2^4%/»"))*tf( 

m'hrnodrn  ktUjSkjiVj  n-0 

by  using  (2.2).  It  follows  that 


(o,-u,4  I)  •••  (»,-«, 4 1) 


V(j)) 


| -v(j>|  s i E | Z 

m h mod  m ’fcj  UjSkjSVj  n*0 

+ I <g«"M«4l)  "•  <V»-U'+l)  . V(J)| 

* Bl  * 

For  fixed  4 - (A|.  • • • . kt)tZs  we  have 

I 2 «<4‘4/«)|  ■ I 2 *(4*4/w)|  » n I e(Aj*ym)|  . 

'kiUjikjZVj  'k}OSkjSVj-Uj  Jm  1 A^-0  1 


VCJ  ' m 

2 ^jkjlm)  ■ Vj-Uj+  I £ m ■ -rr- 


if  hj  m 0 mod  m,  and 


'Aj  - o 


KA ,,  m) 
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. ,_J  \eih/vj-uJ+lVm)-i\  ^ 2 I 

I*  . Mjkjlm) | - |*(6 Jm)-  l]  - |*<6/m)-l|  " n||Aym||  " rfKJTm) 


m 


If  hj  e 0 mod  m,  and  so 


.,„/v v/4'4/")|  s 0 m Hta- 


(2.5) 


In  order  to  estimate  the  second  term  on  the  right-hand  side  of  (2.4),  one  shows  first  by 
induction  on  s that 


IV'W  AI  * £ JVM 

jm  i 

whenever  0<7j,  bji  I for  1 */<  j.  Consequently, 


(w,-u,+  1)  •••  (w,-u,+  I) 


m 


-V(J) 


n^1  ' 


S 

sS 


>1 


m 


-n 


}•  i 


Vj-Uj*  1 


m 


-(0 j-a j)  . 


From  the  definition  of  Uj  and  Vj  it  follows  that 


^ UJ  . 1 . i s 1 ^ VJ  ^ o 

aJS-S<a^m  >nd  BS  mim<6l 


so  that 


v,-u,+  \ „ i 

- 0raj>  < £ for  I <J£s. 


m 


Therefore. 


(i>,-u,  + l)  •••  (t>,-u,+  I) 


m 


-V<J>I  < S’ 


and  by  combining  this  with  (2.4)  and  (2.5),  we  arrive  at 

* 


N-\ 

I I I 

h mod  m ^ n-0 

In  view  of  (2.3),  this  Inequality  holds  for  all  J,  and  by  forming  the  supremum  over  J on  the 


left-hand  side,  we  obtain  the  desired  inequality  for  I 


The  latter  part  of  this  proof  Indicates  that  it  Is  somewhat  unfair  to  use  arbitrary  intervals  J 
when  computing  the  discrepancy.  The  points  <yn>  lie  on  the  coordinate  lattice,  and  we  may 
as  well  assume  that  J is  an  Interval  of  the  type  [a,/m,  6,/m)  x . . . x [ajm,  bjm).  This 
assumption  disposes  of  the  term  s/m  In  the  statement  of  Lemma  2.1. 

11 


Discrepancy  and  Exponential  Sums 


Lemma  2.2.  For  any  integer  mi 2 we  have 

I 


2 -tt — t <(-logm*^) 

m^  n b) 

Proof.  [Niederreiter78]  Since 

2 

h mod  m 


'(S’ 


Ki.  m) 


r{h,  m) 


h mod  m 

it  suffices  to  estimate  the  sum  on  the  right-hand  side.  We  have 

* 


)*. 


2 

h mod  771 


fl  * ^ i n Cm/2j  . 

r-U-  y CSC  nil  — II  £ U - y.  CSC  — , 
r\ft,  m)  r,  _ m " m m ’ 


' h mod  m 

and  by  comparing  sums  with  integrals  we  get 


[m/2]  [m/2]  . rlm/2] 

V*  nn  n ^ nh  ^ n | nx  , 

2j  CK  --  - CSC  - ♦ 2j  esc  — < CSC  - ♦ I CSC  — dx 

h-i  m m h*2  m rn  Ji  m 

Xn/2 

. ..  n n, _ n ^ n m,  2m 

. CSC  t dt  m CSC  — ♦ — k)g  cot  jr—  < CSC  — ♦ — lOK  . 

fm  Tfi  ft  citi  in  ft  n 

For  mi 6 we  have  (min)  sin  (n lm)i  (6/n)  sin  (n/6),  hence  sin  (nlm)i$lm.  This  Implies 


tm/a] 


2 **  ^ m*(  5 - ^ *°g  5)”  for  m2 6, 


I I 


h m 1 


and  so 


[m/2] 


V n/»  . m. 
Z#  CJC  — < —log 
",  m n 6 


m* 


for  m2 6. 


h - 1 


This  last  inequality  is  easily  checked  Tor  m - 3,  4.  and  5,  so  that 


h mod  m 

For  m ■ 2,  Lemma  2.2  is  shown  by  inspection.  I 


The  papers  lNiederrelter72J,  [Niederreiter74],  and  [Niederrelter76]  contain  further  theory 
of  discrepancy  and  exponential  sums. 
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CHAPTER  S 

Reduction  to  a Weaker  Criterion 


Chapter  I introduced  the  linear  recurrences  we  are  investigating  and  gave  some  of  the 
known  results  on  the  length  of  the  period.  Chapter  2 defined  discrepancy  as  a measure  of 
the  distribution  of  a set  of  points  and  presented  Niederreiter’s  lemma,  which  bounds  the 
discrepancy  by  an  exponential  sum.  Computing  the  exponential  sum  is  straightforward  but 
costly;  even  for  small  cases  the  amount  of  computation  becomes  prohibitive.  Since  we  are 
interested  in  the  exponential  sum  mainly  as  a bound  on  the  discrepancy  and  hence  as  an 
indicator  of  the  goodness  of  the  distribution  of  the  sequence,  in  this  Chapter  S we  bound 
the  exponential  sum  by  a function  involving  the  number  of  zeroes  occurring  in  a related 
sequence.  The  number  of  zeroes  can  be  much  larger  than  expected  and  the  discrepancy  of 
the  original  sequence  will  still  approach  zero.  The  question  of  equidistribution  of  the 
original  sequence  is  thus  reduced  to  a weaker  distribution  criterion. 

In  the  case  of  higher  order  congruences  we  cannot  expect  the  discrepancy  to  be  less  than 
m's,  since  some  points  are  never  generated  by  the  recurrence,  for  example,  no  integer 
congruent  to  6 mod  8 will  ever  appear  in  the  Fibonacci  sequence.  The  problem  in  proving 
small  discrepancy  lies  in  showing  that  the  values  which  occur  more  often  than  expected  do 
not  occur  too  often.  Intuition  for  this  problem  comes  from  considering  recurrences  based  on 
primitive  trinomials  with  unit  coefficients  and  looking  at  the  carries  that  occur  when  the 
addition  in  the  recurrence  is  performed  in  radix-2  positional  notation.  The  distribution  of 
carries  should  say  something  about  the  distribution  of  the  sequence. 

In  order  to  study  the  relationship  between  carries  and  distribution  of  digits,  let  us  consider 
the  top  S bits  and  the  bottom  S bits  of  <Fn  mod  64>.  The  period  of  <Fn  mod  8>  is  12  and 
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the  period  of  <Fn  mod  64>  is  96  In  the  matrix  formulation  of  Eqs.  (1.3)  and  (1.4)  we  have 

A - (°  J)  and  A12  - 233)  “ (q  + 8(is  29^  “ , + 8B-  Tak,nS  the  difference  of  two 

vectors  12  apart,  we  see  that  jtM2- Jt  - A - (A ,2- 1)^  - 8Bjr  and  in  general 

” fcr  ■ 8jBjt  mod  64  . (3.1) 


Breaking  the  period  of  <F  n mod  64>  into  8 blocks  of  12  and  considering  each  block  as  a 
point  in  12-space,  the  points  are  a+8tb  mod  64,  where  a,  t « Z12  and  OitiT.  The 
coordinates  of  g are  the  first  12  terms  of  <Fn  mod  64>.  The  coordinates  of  Jj  satisfy  the 
Fibonacci  recurrence  bn  - 6n_, +6n_2  mod  8,  since  they  are  the  difference  of  two  Fibonacci 
sequences.  Here  are  two  tables  illustrating  the  relationship  implied  by  Eq.  (3.1).  The  first 
gives  <Fn  mod  64>  and  the  second  gives  <Fn  mod  64>  where  <Fn>  satisfies  the  Fibonacci 
recurrence,  but  with  initial  conditions  F0-2,  F,«5.  The  entries  In  each  table  are  octal 
integers. 


Table  1.  Octal  values  of  F,2J+t  mod  64 


j 

0 

t 0 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

II 

00 

01 

01 

02 

03 

05 

10 

15 

25 

42 

67 

31 

1 

20 

51 

71 

42 

33 

75 

30 

25 

55 

02 

57 

61 

2 

40 

21 

61 

02 

63 

65 

50 

35 

05 

42 

47 

II 

3 

60 

71 

51 

42 

13 

55 

70 

45 

35 

02 

37 

41 

4 

00 

41 

41 

02 

43 

45 

10 

55 

65 

42 

27 

71 

5 

20 

11 

31 

42 

73 

35 

30 

65 

15 

02 

17 

21 

6 

40 

61 

21 

02 

23 

25 

50 

75 

45 

42 

07 

51 

7 

60 

31 

II 

42 

53 

15 

70 

05 

75 

02 

77 

01 

20 

50 

70 

40 

30 

70 

20 

10 

30 

40 

70 

SO 
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Table  2. 

Octal  values  of 

F \zj*t  moc*  ®4 

j 

0 

t 0 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

II 

02 

05 

07 

14 

23 

37 

62 

21 

03 

24 

27 

53 

1 

02 

55 

57 

34 

13 

47 

62 

31 

13 

44 

57 

23 

2 

02 

25 

27 

54 

03 

57 

62 

41 

23 

64 

07 

73 

3 

02 

75 

77 

74 

73 

67 

62 

51 

S3 

04 

37 

43 

4 

02 

45 

47 

14 

63 

77 

62 

61 

43 

24 

67 

13 

5 

02 

15 

17 

34 

53 

07 

62 

71 

53 

44 

17 

63 

6 

02 

65 

67 

54 

43 

17 

62 

01 

63 

64 

47 

33 

7 

02 

35 

37 

74 

33 

27 

62 

II 

73 

04 

77 

03 

00 

50 

50 

20 

70 

10 

00 

10 

10 

20 

30 

50  row  difference 

This  observation  can  be  applied  directly  to  the  evaluation  of  the  sum  in  Niederreiter’s 
Lemma.  Let  m - p2®,  <j„>-  successive  tuples  of  <yn  mod  m»,  and  N « length  of  one 
period  mod  m.  Then 


Sifajp**) 


11*0 


W/p®-t  pP-i  W/p“-l  p®-1 

£ £ ‘(H<ij+pattypZa)-  £ *&'iJipza)'Z  ^ypa)(3.2) 

J-0  t-0  J-0  f-0 


The  inner  sum  is  a geometric  series  with  ratio  Ki’fe/P®)-  In  fact,  the  sum  is  zero  if  **  0 
mod  pa,  and  is  pa  if  • 0 mod  pa.  Also,  the  sum  depends  only  on  4 mod  p®,  which  is 
4 mod  n,/2. 


Next  we  observe  that  <4’J„>  satisfies  the  same  recurrence  as  <yn>.  To  show  this,  let 
Jn-blrr  Then 

A-J  A-1  A-J 

2 aJf'n*J  " 2 afi'ln>J  -4*2  a£n*J  m * 7n*A  (3.3) 

J»  0 J*0  7*0 

Furthermore,  if  jn  contains  an  element  relatively  prime  to  p,  and  If  ssk,  the  g.c.d.  of  <4*Jn> 
Is  equal  to  the  g.c.d.  of  4;  that  is,  <4*J„>  does  not  degenerate  into  a sequence  in  which  all 
elements  are  divisible  by  a power  of  p higher  than  the  highest  power  of  p dividing  4.  It 
suffices  to  prove  the  case  for  which  the  g.c.d.  of  4 I*  !•  Extend  4 to  a 4-tuple  (Ap  . . . , 
,)  by  adding  k-s  zero  coordinates,  and  define  the  polynomial  H(x)  by 
H(x ) - AqX0*  • • • + ^a-i*A  * • Then  the  polynomial  associated  with  <&•%„>  by  Eq.  (1.2)  is 
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H(x)Y(x)  (mocid  p,  fix))  . In  the  field  Zp[xV{f{x))  both  H(x)  and  K(x)  are  nonzero,  hence 
tiirji  |Modnct  is  also  nonzero.  Therefore  p does  not  divide  all  of  the  coefficients  of 
fHx)Y(x)  ami  thus  <h-yn>  contains  an  element  relatively  prime  to  p.  The  foregoing 
observations  about  <b-yn>  allow  us  to  get  rid  of  vector  operations  and  return  to  plain 

ii« 

integer  sequences.  Let  aj  * h-dj,  bj  « h-bj,  yn  * h>yn.  Then 

N- 1 N- 1 N/pn- 1 prt- 1 

£ * £ 'tfjp2*)  - £ *wjpa)  <3/« 

n- 0 n- 0 J*  0 t»0 

I iie  full  sum.  and  hence  an  estimate  of  the  discrepancy  for  sequences  which  are  2a  bits 
wide,  depends  on  the  occurrences  of  0 in  <bj>,  which  is  only  a bits  wide.  Again  we  note 
that  the  inner  sum  is  a geometric  series  with  a term  ratio  of  t(b Jpa),  and  its  sum  is  0 if 
bj  r (i  mod  pa  or  pa  if  bj  z 0 mod  pa.  Therefore 

N/Pa- 1 prt-l  N/pa- « 

£ e(a'jp?a)Z  eitbjip")*  p<*  £ 4*Jp**)  (3.5) 

J-0  f0  J*  0 

bjso  mod  pP 


Taking  the  absolute  value, 


I 

i pa  £ «<«>ao)i 

J'O 

bj  = 0 mod  pn 


/V/p?-l 

s/>a  £ 


bj  = 0 mod  pa 


pa'(number  of  zeroes  in  <bj  mod  pa>). 


'I  lie  number  of  zeroes  in  <bj  mod  pa>  will  be  large  if  a high  power  of  p divides  h,  since 
bj  *-  d’bj-  To  analyze  this  situation,  let  us  go  back  to  the  original  exponential  sum  of 
Niederreiter’s  lemma  and  choose  m - pa.  Set 

* T<pfl,)-l 

*«> • 2 — r ~ E <<‘-j wn 

> T(f‘  ) 

where  T (pa)  is  the  period  modulo  pa,  and  let  a 0 be  an  exponent  for  which 
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T (f>ao*J)  pJr{pao)  for  all^iO.  Then  for  cr>cr0 

or  T (pa)-i 

2 2 ~ ~~  2 f<!L‘XrJPa) 

Pmodp"  '<M">  *<*“> 

gcd(h)  «■  pa~d 

a T(pa)-i 

’ £ » ii  »r  w")  «**>  »?° 'w°  "'I” ' ('l>“ d>> 

gcd(h)  « pa‘d 


5(a0-.)+  £ r(<,-d>  • 

rf'-Or0  P mod  pd  r(-'  P > 

gcd(h)  « t 


Thus  for  a >o0 

T(pn)-i 

•S(n)  = p~5S(a- 1)+  2 2 Pa)  (3  8) 

hmodpafVL.pa ) npa)  - o 

gcd(h).  I 

Collapsing;  the  chain  of  Inequalities  back  to  the  statement  of  Niederreiter’s  lemma,  for 


£ <«•>„/ 

r r , T (,*>)  ft  -31" 

*ed(P) « I 


< JL  ♦ ♦ V ! £l_ 

^ pmodp^P**) 

gcd(h)  - 1 


— • — Knumber  of  zeroes  In  mod  pa>) . (3  9) 
2«\  w/,2as  * r 


For  or>o0  the  g.c.d.  of  the  coni|K>nents  of  bj  Is  I,  by  Lemma  1.4.  Thus  by  an  earlier 
remark  the  g.c.d.  of  the  <bj>  is  I,  because  b'j  - h-bj  and  the  sum  (3.9)  is  restricted  to  those  h 
fur  which  gc.d.(A)  - I.  Replacing  the  number  of  zeroes  in  <bj  mod  pa>  by  the  maximum 
number  of  zeroes  in  a nondegcncrate  cycle  mod  pa  and  applying  the  bound  of  Lemma  2.2 
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gives 

„ i 5(2or- 1) 

T(/>  ' pZCt  p* 


pCt 

(f.2a)5*  - ' •(maximum  number  of  zeroes  in  nondegenerate  cycle  mod  pa) . (3.10) 

T </>*») 


for  some  c , which  depends  on  p but  not  on  a Remembering  that  T (pZct)  is  proportional  to 
p?ct,  we  see  that  the  discrepancy  in  s<k  dimensions  will  approach  zero  if  the  maximum 
number  of  zeioes  in  a cycle  modulo  pa  is  o(ct~spa).  This  is  a very  weak  condition.  If  the 
elements  of  the  cycle  are  evenly  distributed  then  the  expected  number  of  zeroes  is  a 
constant.  Equation  3 10  shows  that  the  discrepancy  will  tend  to  zero  even  if  the  number  of 
zeroes  in  a cycle  Is  exponentially  increasing,  as  long  as  the  rate  of  increase  is  />-t  for  some 

€ N 0. 
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CHAPTER  4 


Analysis  of  the  Number  of  Zeroes  in  a Cycle 


The  result  of  Chapter  3 expresses  a bound  on  the  discrepancy  in  terms  of  the  number  of 
zeroes  occurring  in  sequences  satisfying  the  recurrence.  (See  [Hall38b]  for  an  early 
application  of  a similar  bound.)  In  this  chapter  we  will  investigate  the  number  of  zeroes 
appearing  in  any  cycle,  and  try  to  bound  it  from  above  by  a function  which  is  small 
enough  to  force  the  discrepancy  to  zero  as  the  modulus  increases. 

For  the  Fibonacci  sequence  modulo  2°  the  bound  is  particularly  small. 

Theorem  4.1.  At  most  two  zeroes  appear  in  a cycle  satisfying  the  Fibonacci  recurrence 
yn, 2 = >„♦!  +y„  modulo  2®  when  y0  and  y,  are  relatively  prime. 

Proof.  If  a zero  appears  at  all  then  shift  the  cycle  so  that  the  zero  is  the  first  element. 
Then  y0  s 0,  y,  e a,  and  a must  be  odd  since  y0  and  are  relatively  prime.  Thus  the 
cycle  is  merely  a multiple  of  the  Fibonacci  sequence,  where  y0  - 0,  y,  - I.  Modulo  8 the 
Fibonacci  sequence  is  01  123505527  1.  By  inspection  there  is  one  zero  mod  2,  one 
zero  mod  4,  and  there  are  two  zeroes  mod  8.  We  will  show  by  induction  that 
F^^a-z  - 2 aq  and  F 3.20-2.,  - 1 + 2°*^  where  or *3  and  q,  r are  odd  integers.  This  is 
true  by  inspection  for  o - 3.  Assume  that  it  is  true  for  a - j.  Then  using  the  relations 
-F2„+F2n,t  and  F2n  - Fn(2Fn,  ,-F,)  we  have  - ** '*1  + *Mr-2'-; lq) 

and  F^2j-i  4 , - I + 2^r  + 22J*2r2  + 2%2,  which  is  the  property  for  a - j+  1.  Therefore  the 
period  doubles  and  the  number  of  zeroes  remains  constant  as  the  modulus  doubles.  I 

We  can  use  this  fact  to  estimate  the  discrepancy  in  two  dimensions  of  a complete  cycle  of 
the  Fibonacci  sequence  modulo  2®.  Assume  that  S(a)<caz 2*a  2 in  the  notation  of  Eq. 
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(3.8).  For  ot>€  we  have  7/5  + (2/n)log  2°  <log  2a.  By  Eq.  (3.9),  S(a+  l)s 

2'25(a)+2*21/2(a+ l)22‘<a+,)/2  - ca22'a/2‘2  + (a+l)223/2  ' <a+1)/2 

c(a  + \ f2'(a*i)/z  ( 2‘3/2(-^-r)2  + c'123/2  ) < c(a+  |)22’(a+,)/2  if  c>8/(23/2-  I).  If  we 

a+l 

also  select  c large  enough  so  that  5(6)  satisfies  the  inequality  then  a simple  induction  will 
establish  the  bound  for  ai 6 Therefore  Dws2,‘a  + ca22'a/2. 

For  recurrences  of  higher  degree  than  the  Fibonacci  sequence,  the  number  of  zeroes  in  a 
cycle  can  increase  as  the  modulus  increases. 

Theorem  4.2.  If  * > 2 and  aia0  then  the  maximum  number  of  zeroes  In  a nondegenerate 
cycle  modulo  pZa  is  at  least  pa. 

Proof.  Construct  k initial  elements  yv  ...,yk.,  a s follows.  Set  y0sO  mod  p2*3.  Choose 

yx yk, j mod  pa  so  that  *,Bj!0  s 0 mod  pa,  where  B Is  the  matrix  defined  in  Eq.  (1.4) 

and  gj  - (I,  0,  0,  ....  0)  is  a unit  vector  with  a I in  the  first  coordinate.  Since  y0  has  been 
specified,  the  constraint  on  y,,  ....  yk.x  is  one  equation  in  * - 1 unknowns.  For  * - 2 (as 
in  the  Fibonacci  sequence),  yx  would  have  to  be  0 mod  pa  and  g.c.d.(y0,  yx)  would  be  pa. 

For  *>2,  well-known  methods  guarantee  the  existence  of  a nontrivial  solution  for  yx 

yk.y  It  is  easy  to  see  that  y n s g,A/rj0  s £x(l+jpaB)y0  s g,.j0+  ty)paBy0  s 0 mod 
p 20  for  0<j<pa  and  that  y„  are  distinct  elements  of  a cycle  modulo  p Za.  I 

In  order  to  obtain  a good  bound  on  the  discrepancy,  it  is  necessary  to  show  that,  as  the 
modulus  goes  from  pa  to  pa*x,  the  number  of  zeroes  increases  by  a factor  which  is 
eventually  less  than  p.  We  will  try  to  show  that  the  number  of  zeroes  eventually  increases 
at  a rate  no  faster  than  p(k‘z'),(K~x\  where  k is  the  degree  of  the  recurrence. 

The  fundamental  idea  will  be  counting  in  two  different  ways  the  total  number  of  zeroes  in 
the  cycles  corresponding  to  all  of  the  possible  initial  conditions.  Consider  a large  tabular 
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array  in  which  each  cycle  appears  as  a row.  The  first  column  of  the  array  contains  the  first 

element  of  each  cycle;  the  second  column  is  composed  of  the  second  element  from  each  cycle, 

and  so  on.  Let  b,  be  the  number  of  zeroes  in  the  fh  row  of  the  array.  Then  2 ( /) 

J 

counts  the  number  of  (-tuples  of  zeroes  in  the  array,  with  all  elements  of  a tuple  required  to 
appear  in  the  same  row.  The  number  of  (-tuples  of  zeroes  can  be  counted  another  way. 
For  each  (-tuple  of  column  indices,  count  the  number  of  rows  which  have  zeroes  in  all  of 
those  ( columns.  Since  the  two  methods  of  counting  must  agree,  a bound  achieved  by 
considering  one  counting  method  can  be  applied  to  the  other  counting  method. 

To  count  the  zeroes  appearing  in  all  of  the  cycles,  we  must  first  know  how  many  cycles 
there  are.  Each  cycle  is  determined  by  its  first  k elements  and  each  element  can  take  on  pa 
values,  so  there  are  pak  possible  cycles.  Some  of  these  cycles  are  isomorphic  under  cyclic 
shift,  and  of  course  the  number  of  zeroes  in  cycles  which  are  cyclic  shifts  of  one  another  is 
the  same.  To  reduce  our  counting  effort  by  applying  knowledge  of  the  isomorphism,  we 
should  divide  by  the  period.  A cycle  in  which  all  the  elements  are  divisible  by  p has  a 
shorter  period  than  a cycle  containing  an  element  relatively  prime  to  p,  so  it  is  necessary  to 
count  the  cycles  according  to  the  highest  power  of  p dividing  all  elements.  Let  aj  be  the 
number  of  cycles  for  which  p * is  this  highest  power.  Then  by  counting  all  the  cycles  we 
have  a0+a,  + " ' +aa  " P°K-  ^ P°  is  **^8*  enw,gh  *°  that  the  period  is  multiplied  by  p 
as  the  modulus  increases  from  pa  to  pa  + then  by  considering  what  happens  when  the 
zero  unit’s  digit  is  removed  from  a cycle  in  which  all  the  elements  are  divisible  by  p,  we  see 
that  <»,+  •••  +oa  - p(a'l)k.  Thus  a0  - pak-p<a'i)k  . (pk-\)p(a't)k.  Considering 
only  those  cycles  containing  an  element  relatively  prime  to  p and  allowing  for  cyclic  shifts, 
there  are  (pk - \)p(a'i)kKpk - l)pa','(e',)  - cycles. 

Some  of  these  cycles  are  isomorphic  under  multiplication  by  a number  prime  to  p.  Cycles 
which  are  multiples  of  each  other  have  the  same  number  of  zeroes.  In  fact,  a cycle  may  be 
non-trlvlally  Isomorphic  to  Itself  under  such  a multiplication.  For  instance,  in  the 
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Fibonacci  sequence  modulo  8 the  second  half  cycle  0 5 5 2 7 1 is  5 times  the  first  half 
cycle  0 11  2 3 5.  In  this  situation  the  number  5 is  called  a multiplier.  If  a cycle  has  a 
multiplier  then  obviously  the  number  of  zeroes  in  the  cycle  is  a multiple  of  the  number  of 
zeroes  appearing  in  the  first  partial  cycle,  which  is  called  a block  [Ward33l  Multipliers 
complicate  some  of  the  later  analyses.  Any  constant  number  of  multipliers  can  be  tolerated, 
but  it  is  preferable  not  to  deal  with  multipliers  at  all.  Therefore  we  must  be  able  to  detect 
when  a cycle  has  a multiplier. 

Suppose  that  a is  a multiplier  after  n iterations  of  the  recurrence,  i.e.,  yn  s ay0,  yn+ , = ay ,, 

....  Let  g(x)  be  the  polynomial  corresponding  to  the  initial  conditions  (y0,  y, yA,). 

Then 

xng(x)  = ag(x)  (modd  pa,flx)) 

(xn-a)g(x)  = 0 (modd  pa,fix )) 

(xn  - a)g(x)  - f[x)u(x)  + pav(x)  for  some  u(x),  v(x) . 

Since  fix)  is  primitive  in  Zp(x],  either /(x)  divides  xn-a  or  /(x)  and  x"-o  are  relatively 
prime.  Suppose  they  are  relatively  prime.  Then  by  Hensefs  Lemma  we  can  find  e(x),  d(x) 
such  that  (x"-a)c(x)+/(x)rf(x)  s I mod  pa.  Multiplying  by  g<x)  we  find 
(xn  - a)g(x)c(x)  + f{x)g(x)d(x)  s g(x)  mod  pa 
/(x)u(x)c(x)+  pav(x)c(x)+ f(x)g(x)<{(x)  s g(x)  mod  pa 

0 = gix)  (modd  pa,flx)). 

Thus  a would  be  a multiplier  for  the  trivial  cycle  (0)  only.  Consider  the  other  possibility, 
that  J[x)  divides  x"-a  in  Zplx]  so  that  x"  s a (modd  p,f,x)).  The  multiplier  a must  be 
relatively  prime  to  p or  else  the  recurrence  could  be  run  backwards  to  deduce  that  all  of  the 
initial  values  were  non-prime  to  p.  Let  j be  the  order  of  a mod  p.  Then  x*n  s aJ  e I 
(modd  p,  f{x)).  Thus  Tip)  divides  Jn.  Noting  that  p-  I divides  pk-  I and  that  j divides 
p-  I,  we  have  Tip)lip-\)  divides  n.  Specializing  to  p - 2,  T(2)  divides  n.  If  we  consider 
the  case  in  which  n is  the  smallest  number  of  iterations  producing  a multiplier  then  n 
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divides  T (pa).  For  p - 2 these  two  divisibility  conditions  restrict  n to  the  form  2a'c(2*-  1), 
and  thus  a is  a 2J  root  of  unity  for  some  J.  The  powers  az,  a4,  a8,  ...  must  also  be 
multipliers,  and  in  particular  one  of  the  square  roots  of  I must  be  a multiplier.  The  square 
roots  of  i mod  2°  are  +i,  -i,  2a'1+l,  and  2a‘,-l.  The  case  of  +1  is  trivial.  If  -I  or 
2a‘1-i  were  a multiplier,  it  would  correspond  to  2'1T(2a)  - T(2a*1)  iterations, 
contradicting  the  known  period  mod  2a~1.  Therefore  2a‘,+  l must  be.  a multiplier 
corresponding  to  one-haif  the  period.  In  matrix  formulation  of  Eq.  (1.3)  and  (1.4), 

AT/2x  e (2a"1+  l)x  mod  2° 

(2°",B  + I)*  e (201*1  + l)x  mod  2° 

2a'1(B-I)x  e Omod  2° 

Thus  B-I  is  singular  mod  2.  For  recurrences  of  moderate  degree  this  is  not  hard  to 
determine.  In  the  nonsingular  cases  there  can  be  no  multipliers. 

Back  to  considering  cycles  isomorphic  under  multiplication  by  a number  relatively  prime  to 
p,  we  now  assume  that  the  recurrence  itself  has  no  multipliers.  We  can  divide  the  number 
of  cycles  by  $(pa),  which  in  the  case  p ■ 2 leaves  c,2*a",)(A‘,)/2a'1  - c,2(a',^*‘2*  cycles 
containing  an  odd  element  and  non-isomorphic  under  cyclic  shift  and  multiplication  by  a 
constant. 

Knowing  the  total  number  of  cycles,  we  begin  counting  the  xeroes  by  counting  the  rows 
which  have  zeroes  in  specified  columns. 

Thtorm  4.3  [Hal138aX  There  exists  a sequence  <yn>  satisfying  the  recurrence  (1.1)  and  not 
Identically  zero  modulo  pa  for  which  yn  * 0 mod  pa  for  A - 1 arbitrary  values  of  n. 

Proof.  Let  the  arbitrary  values  of  n be  n,,  ...,  nk.v  Write 

<y„>  - <e0w„+ciwn.i+  •••  +ek.i«yn*a-i>  where  **1?  u the  un,t  *®9U«Ke  with  initial 
conditions  (0,  0,  ...,  0,  I)  and  the  c’s  are  to  be  determined  by  the  k-l  congruences 
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^yn i s yn  = •••  s yn * ^ mod  pa  . These  are  A-  I homogeneous  linear  congruences  in  the 
A variables  c & . . . , ck,  so  there  must  exist  a solution  in  which  not  all  the  c's  vanish  and  yn 
does  not  vanish  identically.  I 

Thus  at  least  one  cycle  has  zeroes  for  each  set  of  A~l  arbitrary  positions.  To  calculate  the 
exact  number  of  cycles  having  zeroes  in  the  specified  positions,  we  need  to  know  the  rank 
of  the  coefficient  matrix.  Suppose  that  the  rank  was  A-l.  Then  there  would  be  one  free 
parameter  in  the  solution  set  and  the  number  of  solutions  would  be  proportional  to  pa  as  a 
varied.  Varying  the  parameter  would  merely  generate  a solution  which  was  a multiple  of 
the  previous  solution,  so  that  the  corresponding  cycles  would  also  be  multiples  of  each 
other.  Thus  if  the  rank  of  the  coefficient  matrix  was  A-  I then  there  would  be  a constant 
number  of  cycles  having  zeroes  in  the  designated  columns.  To  take  advantage  of  the  cyclic 
shift  isomorphism,  we  can  demand  that  first  column  be  one  of  the  columns  containing  a 
zero.  This  leaves  A -2  other  columns  which  can  be  specified,  so  altogether  the  number  of 
systems  of  simultaneous  homogeneous  linear  equations  we  are  considering  is  proportional  to 
pa(.k-Z)  We  f,ave  Jeen  that  jf  t^e  coeffjdjpt  matrices  are  of  rank  A - I then  each  system 
corresponds  to  a constant  number  of  cycles.  Letting  bj  be  the  number  of  zeroes  in  the  jth 
row,  this  argument  shows  that  £ where  the  sum  is  over  nondegenerate 

cycles  with  a zero  in  the  first  position.  This  Implies  that  A**1  iczpa^k'2^  for  each  J,  and 
hence  the  maximum  number  of  zeroes  in  a cycle  mod  pa  would  be  bounded  by 
c3pa('k'2'^,<'k for  some  constants  c#  c 3.  This  bound  is  good  enough  to  force  the 
discrepancy  to  zero  as  a increases. 

For  recurrences  of  degree  3 that  have  no  multipliers,  the  coefficient  matrix  is  of  rank 
2 - A - I.  The  rank  cannot  be  zero  because  no  A consecutive  terms  vanish,  and  if  the  rank 
were  I then  the  second  row  would  be  a multiple  of  the  first.  Hence  the  preceding  argument 
proves  the  following  theorem. 
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Theorem  4.4.  If  J{x)  is  a primitive  polynomial  of  degree  3 in  ZJlx)  then  the  discrepancy 
over  an  entire  period  of  f-tuples  generated  by  (l.l)  tends  to  zero  as  a tends  to  infinity,  for 
tii. 

The  cases  in  which  the  coefficient  matrix  determined  by  the  k - 1 column  positions  is  of 
rank  A - 1 - 1 for  t> 0 generate  a number  of  essentially  distinct  cycles  proportional  to  pat. 
Since  the  matrix  is  not  of  full  rank  there  is  some  non-trivial  linear  combination  of  the  rows 
which  is  the  zero  vector.  Using  the  isomorphism  (1.2)  between  sequences  and  polynomials, 
the  same  non-trivial  linear  combination  of  the  powers  of  x corresponding  to  the  designated 
column  positions  must  be  the  zero  polynomial  modd  pa,J{x).  This  means  that  we  are  now 
interested  in  the  number  of  polynomials  with  k-t  terms  (one  of  which  is  the  constant  term 
c0x°)  which  are  congruent  to  zero  modd  pa,  fix).  We  can  choose  k-t-  l exponents  of  x 
from  a number  of  possible  exponents  proportional  to  pa.  The  coefficient  for  each  of  the 
k - 1 terms  can  be  chosen  in  pa  ways,  but  this  gives  pa  times  too  many  choices  because  of 
constant  multiples.  Thus  the  number  of  polynomials  to  consider  is  proportional  to 
pa(k-t-i)pd(k-t)p-a  m pZa(k-t-i)  jj.  t^ese  polynomials  were  evenly  distributed  among  the 
pak  residue  classes  modd  pa,f{x)  then  the  probable  number  of  polynomials  congruent  to 
zero  would  be  pQ(k~2t~2\  an(j  the  number  of  (A-l)  tuples  of  zeroes  would  be  proportional 
to  pa***t'2*.  Summing  for  t from  0 to  A -2  gives  a number  proportional  to  pa^k'2^  and 
therefore  bJ<cApa{k~z)l{k‘'). 

Another  way  to  attack  the  problem  of  bounding  the  maximum  number  of  zeroes  in  a cycle 
modulo  pa  uses  the  matrix  formulation  of  the  recurrence.  Let  A be  the  matrix  (1.3). 
Assume  that  a is  large  enough  so  that  the  period  is  multiplied  by  p as  or  increases  by  I, 
and  that  yn  m 0 mod  pa.  Let  T ■ T (pa)  be  the  period  mod  pa,  let  j - ....  ?„♦*.  j)T 

be  the  A-tuple  of  residues  at  y#  let  B be  the  matrix  (1.4)  so  that  A T - I + p0B,  and  let 
l *JI4*.1)T  - By.  Let  ^°4,Hyfl  and  P^JUt^  If  l<]  then  modulo  pa  *J  * 1 there 

will  be  no  zeroes  with  indices  congruent  to  n modulo  T (pa).  If  UJ  then  as  a increases 
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there  will  eventually  be  ( p - I )pJ  zeroes  in  each  period  with  indices  congruent  to  n modulo 
T (pa).  This  shows  that,  given  any  initial  conditions  over  the  integers,  the  number  of  zeroes 
in  a cycle  generated  from  those  initial  conditions  will  eventually  be  a constant.  The 
problem  is  that  as  the  modulus  increases  the  sets  of  initial  conditions  which  are  allowed 
also  varies;  the  number  of  zeroes  in  a cycle  converges  pointwise  (for  each  set  of  initial 
values  over  the  integers),  but  the  question  of  uniform  behavior  is  as  yet  unanswered. 


Establishing  a nontrivial  bound  on  the  number  of  zeroes  in  a period  of  a recurrence 
satisfying  (l.l)  is  a worthwhile  future  project.  Any  bound  which  is  o (a~kpa)  can  be  used 
in  Eq.  (3.10)  to  show  that  the  discrepancy  of  a full  period  of  A-tuples  tends  to  zero  as  a 
tends  to  infinity.  The  conjecture  below  might  possibly  be  proved  by  establishing  suitable 
bounds  on  the  distribution  among  the  residue  classes  modd  pa,/(x ) of  polynomials  with  k 

. 

or.  fewer  terms  and  degree  less  than  T (pa),  or  by  other  means. 


Conjecture  4.1.  The  maximum  number  of  zeroes  in  a nondegenerate  period  of  a recurrence 
generated  by  Eq.  (l.l)  is  less  than  cp0^*'2^*'1)  for  some  c depending  on  p and  k but  not 
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Practical  Considerations 


In  Chapter  4 we  saw  indications  that,  when  considered  over  an  entire  cycle,  the  tuples  of 
residues  probably  become  equidistributed  as  the  recurrence  Is  computed  with  more  and 
more  bits.  In  any  practical  situation  only  a moderately  large,  fixed  number  of  bits  are  used; 
most  current  machines  use  from  24  to  48  bits  to  store  the  fraction  of  a single-precision 
floating-point  number.  Even  so,  it  Is  highly  unlikely  that  an  entire  period  will  ever  be 
used.  Choosing  p - 2,  a - 24,  and  *>50  implies  a period  of  about  270,  which  Is  in  the 
range  of  the  total  number  of  states  that  all  the  computers  In  the  world  have  ever  been  in 
(220  machines  • 26  years  • 226  seconds/year  • 220  states/machine-second).  A truly  random 
sequence  is  locally  non-random  in  some  places;  such  places  should  not  occur  frequently  In 
the  portions  of  a computer-generated  sequence  which  are  likely  to  be  used.  To  be 
recommended  for  general  use,  a random  number  generator  must  be  efficient  and  easy  to 
implement  on  a variety  of  machines.  This  chapter  addresses  these  practical  considerations 
and  presents  the  results  of  comparisons  with  a standard  linear  congruential  generator 

>n*l-a>n+6mod2a- 

A random  number  generator  based  on  recurrence  (i.i)  can  be  efficiently  implemented. 
Only  the  * most  recent  values  of  y need  to  be  remembered,  and  these  can  be  stored  In  an 
array  which  is  accessed  cyclically.  The  number  of  arithmetic  operations  involved  can  be 
minimized  by  choosing  fix)  to  be  a primitive  trinomial  modulo  2.  (Some  primitive 
polynomials  mod  2 have  been  tabulated  in  [Watson62]  and  [ZierlerSSj.)  If  fix)  is  chosen  in 
this  way  then  the  recurrence  is  yn  • *yH.j±yH. * mod  2°,  which  requires  no  multiplication 
or  division,  only  one  addition  or  subtraction,  and  one  reduction  modulo  2°.  Many 
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applications  require  only  the  fractions  These  fractions  can  be  computed  exactly  in 

floating-point  arithmetic  without  intermediate  integer  computation  and  conversion.  This  is 
done  by  dividing  each  term  of  the  recurrence  by  2°,  so  that  the  computation  is  done  in 
floating-point  numbers  modulo  I.  It  is  necessary  to  prevent  any  floating-point  operation 
from  generating  a result  greater  than  1.0,  which  would  involve  a right  shift  of  the  fraction 
part  and  the  possibility  of  losing  one  bit,  making  the  computation  no  longer  exact.  This 
can  be  insured  by  choosing  the  signs  so  that  J{x)  ■ xk-xk'J+  I and  the  recurrence  is 
yn  “ Jn-j~yn-k  mod  Th,s  ,s  allowed  because  +1  e - I mod  2 and  therefore  x* -**'•*+  1 
is  the  same  polynomial  in  Z2  as  xA+x*"^+  I.  If  the  subtraction  jn_y2a  - produces 

a negative  result  then  we  add  1.0  to  bring  the  value  back  into  the  range  0sjtJI/2a<  I. 
When  performing  this  addition  we  must  remember  to  actually  perform  two  additions  of  0.5 
in  order  to  guarantee  that  no  bits  will  be  lost.  (Computers  such  as  the  Control  Data 
Corporation  6000  series  effectively  do  not  use  any  guard  digits  and  hence  one  bit  can  be 
lost  when  preshifting  Jin/2a  to  align  the  radix  point  before  the  addition.) 


To  illustrate  the  fact  that  a random  number  generator  based  on  (l.l)  can  be  implemented 
with  little  machine  dependence,  here  is  a coding  in  FORTRAN  of  a generator  based  on 
x66-x3t  + I mod  2.  Each  time  XRAND  is  called  it  returns  the  next  random  number. 


FUNCTION  XRAN0 

COMMON  /XRAN0X/  I,  J,  X<55) 

DATA  I /55/ 

I - I -1 

IF  (I  .IE.  0)  I - 55 
J - I - 31 

IF  (J  .IE.  01  J - J + 55 
XU)  - XU)  - X(J) 

IF  (XU)  .IT.  0.0)  XU)  - (XU)  4 0.5)  4 0.5 
XRAN0  - XU) 

RETURN 

EN0 


The  machine  dependence  lies  In  the  initialisation  of  the  first  55  values  of  the  array  X.  The 
initial  values  must  be  chosen  so  that  0s XU)  < I,  2aX(I)  is  an  integer,  and  2aX(I)  must 
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be  odd  for  at  least  one  value  of  I.  One  way  to  do  this  is  to  set  X (1 ) - 1/2°  and  X ( 1 ) ■ 0 
for  2 s I s 55.  However,  this  puts  one  of  the  local  non-random  areas  at  the  start  of  the 
generated  sequence.  It  would  be  better  to  run  the  recurrence  a million  times  (say),  write  out 
the  values  of  X(l),  and  have  some  other  routine  place  these  values  into  X before  using 
XRAND.  Alternatively,  an  auxiliary  linear  congruential  generator  can  supply  the  initial 
values;  see  [Brent 73]. 

The  efficiency  of  generator  XRAND  compares  favorably  with  that  of  standard  first-order 
linear  congruential  generators.  On  the  Digital  Equipment  Corporation  PDP-IO,  XRAND  can 
be  encoded  in  12  instructions,  while  a standard  generator  requires  10  instructions.  (The 
number  of  instructions  is  for  a routine  which  returns  a floating-point  value  between  0 and 
I,  is  guaranteed  to  not  cause  any  interrupts,  and  includes  the  normal  subroutine  linkage.) 
Execution  time  per  call  on  a KL-10  processor  Is  approximately  equal  because  XRAND  uses 
no  multiplication.  On  the  IBM  S/370,  comparable  encodings  are  16  instructions  (46  bytes) 
for  XRAND  and  10  instructions  (27  bytes)  for  a first-order  generator. 

Four  generators  were  compared  using  several  tests.  The  generators  were 

RANDU;  y„  - (2,e+3 )yn.i  mod  231,  y0  - I.  This  generator  is  notorious  for  its 
bad  distribution  in  three  dimensions. 

GOODLC:  y„  - 314 1592653yn.,  + 2718281829  mod  236,  j0  - 0.  This  Is  a standard 
"linear  congruential"  generator. 

ADDLC:  yn  - y„.66-yn.24  mod  227,  y0  - • ••  - J63  - °,  J64  - The  tests  began 
with  - y66636.  This  is  an  additive  generator  based  on  the  primitive  polynomial 
x66-x3,+  I mod  2. 

, BESTX:  y„*xn  XOR  t„  where  x„  - (3I4I592653x„_,  ♦ 2718281829)  mod  236, 

x0  ■ 0,  tn  - 314159270*,,.,  mod  (236-  31),  and  x0  ■ I.  The  number  236-  31  is  the 
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largest  prime  less  than  236,  and  314159270  is  a primitive  root.  Here  XOR  is  the 
exdusive-or  bit  operation. 

These  generators  were  compared  under  the  following  tests.  (See  [Knuth69]  sec.  3.3.2.) 

ID:  A one-dimensional  distribution  test  with  the  interval  [0,1)  divided  into  4096 
equal  subintervals,  and  8 hits  expected  in  each  subinterval. 

2D:  A two-dimensional  distribution  test  with  the  interval  [0,1)  x [0,1)  divided  into 
64  x 64  equal  subintervals,  and  8 hits  expected  in  each  subinterval.  The  ordered 
pairs  used  for  the  test  were  non-overlapping,  i.e.,  (jq,  },),  (y2,  y3) 

3D.  A three-dimensional  distribution  test  with  the  interval  [0,1)  x [0,1)  x [0,1) 
divided  into  16  x 16  x 16  equal  subintervals,  and  5 hits  expected  in  each 
subinterval.  The  ordered  triples  used  for  test  were  (y0 , yv  yz),  (y3.  y4,  yB) 

GAP:  A test  which  considers  the  length  of  consecutive  subsequences  yj,  . . . 

yj*r  in  which  Osa  *yjtr<0 < I for  two  fixed  real  numbers  a and  0,  but  the  other 
y’$  do  not.  This  test  was  performed  with  a - 0,  (5  - 0.5,  and  thus  was  a test  of  "runs 
above  the  mean".  Gaps  of  length  0 through  5 and  gaps  of  length  greater  than  5 
were  counted  until  500  gaps  had  been  tabulated. 

MAX  10:  The  maximum  element  of  blocks  of  10  consecutive  values  was  selected 
until  1500  maxima  had  been  chosen.  The  distribution  of  the  maxima  was  tested 
against  the  theoretical  distribution  function  z10  by  the  Kolmogorov-Smirnov  test. 

RUN:  The  length  r of  consecutive  subsequences  for  which 

>yj*rsyj*r*i  was  tabulated  until  500  runs  had  occurred.  New  runs 
were  started  at  yJtrt2.  and  runs  of  length  5 or  more  were  grouped  together  for  the 
analysis. 

PERMUT:  The  order  relations  among  consecutive  blocks  of  4 values  were 


30 


r 


Practical  Considerations 

tabluated,  with  each  of  the  24  possible  orderings  expected  to  occur  150  times. 

Except  for  the  test  MAX  10,  where  the  Kolmogorov-Smirnov  (KS)  test  was  used,  the 
chi-squared  statistic  was  calculated  using  the  appropriate  probabilities.  Exceptions  to  the 
expected  distribution  were  counted  when  the  calculated  statistic  lay  in  the  5*  tail  at  either 
end  of  the  theoretical  distribution.  Thus  a perfectly  random  sequence  would  be  expected  to 
fail  I OX  of  the  tests  in  each  category.  The  chi-squared  values  themselves  were  tested  in 
groups  of  16  by  the  KS  test  (both  KS+  and  KS-  tests)  against  the  hypothesis  that  they 
came  from  a chi-square  distribution.  Exceptions  were  noted  for  the  5t  tails  at  both  ends  of 
the  KS  distribution. 

Each  test  was  repeated  until  a conveniently  large  percentage  of  the  first  S-220  values  of 
each  generator  had  been  tested.  The  following  tables  summarize  the  results. 


Table  3.  Results  of  tests  on  RANDU 


test 

repetitions 

±5X  tails 

KS  tests 

KS  ±5*  tails 

ID 

96 

12 

12 

3 

2D 

48 

5 

6 

1 

3D 

48 

48 

6 

6 

MAX10 

384 

50 

GAP 

192 

20 

24 

4 

RUN 

192 

26 

24 

6 

PERMUT 

192 

II 

24 

0 

Table  4.  Results  of  tests  on  COODLC 


test 1 repetitions  1 

± St  tails 

KS  tests 

KS  ± 5X  tails 

ID 

96 

10 

12 

3 

2D 

48 

5 

6 

1 

3D 

48 

4 

6 

0 

MAXIO 

384 

39 

GAP 

192 

19 

24 

2 

RUN 

192 

29 

24 

9 

PERMUT 

192 

20 

24 

2 
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Table  5.  Results  of  tests  on  ADDLC 


test 

± 51  tails 

KS  tests 

KS  ± 51  tails 

ID 

96 

13 

12 

0 

2D 

48 

6 

6 

0 

3D 

48 

5 

6 

0 

MAX  10 

384 

34 

GAP 

192 

19 

24 

3 

RUN 

192 

29 

24 

8 

PERMUT 

192 

17 

24 

3 

T able  6.  Results  of  tests  on  BESTX 


test 

1 repetitions 

± 51  tails 

KS  tests 

KS  ±51  tails 

ID 

15 

12 

0 

2D 

6 

6 

1 

3D 

3 

6 

0 

MAX  10 

384 

41 

GAP 

192 

16 

24 

3 

RUN 

192 

21 

24 

12 

PERMUT 

192 

13 

24 

0 

The  tests  confirm  the  bad  three-dimensional  distribution  of  the  values  generated  by 
RANDU.  All  generators  had  difficulty  with  the  RUN  test  The  code  for  this  test  was 
carefully  examined  for  systematic  errors,  but  none  were  found.  The  exceptional  chi-square 
values  tended  to  be  extremely  small  (less  than  1.0)  or  Just  above  the  upper  5*  tall  cutoff. 
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Generator  ADDLC  compares  favorably  with  GOODLC  and  BESTX  In  these  tests.  To 
the  extent  that  this  testing  procedure  Is  valid  for  a particular  task  requiring  random 
numbers,  ADDLC  can  be  recommended  as  an  acceptable  generator. 
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